/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Contains declarations of the wd01_grain_model class. */

// $Id$

#ifndef __wd01_grain_model__
#define __wd01_grain_model__

#include <vector>
#include "grain_collection_grain_model.h"
#include "thermal_equilibrium_grain.h"
#include "grain_size.h"
#include "misc.h"

namespace mcrx {
  template <template<class> class, typename> class wd01_grain_model;
  template <template<class> class chromatic_policy, typename T_rng_policy> 
  class wd01_grain_creator_functor;
};


/** This class implements the Weingartner & Draine 01 grain model,
    with size distributions and grain data from that paper. The class
    methods are not reentrant and can not be used by multiple
    threads. For the temperature calculation, you must make a copy of
    the object for each thread.  */
template <template<class> class chromatic_policy, 
	  typename T_rng_policy = mcrx::global_random> 
class mcrx::wd01_grain_model : 
  public mcrx::grain_collection_grain_model<chromatic_policy, T_rng_policy> {

public:
  wd01_grain_model (const Preferences& p, const T_unit_map&);
  wd01_grain_model (const wd01_grain_model& g) :
    // need to make sure array members are copied. We make the copies
    // thread-local, since the methods are not reentrant anyway.
    grain_collection_grain_model<chromatic_policy, T_rng_policy>(g) {};

  virtual boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> > 
  clone() const {
    return boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> >
      (new wd01_grain_model(*this));};
};



template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
mcrx::wd01_grain_model<chromatic_policy, T_rng_policy>::
wd01_grain_model (const Preferences& p,
		  const T_unit_map& units) :
  grain_collection_grain_model<chromatic_policy, T_rng_policy>()
{
  const std::string data_directory =
      word_expand(p.getValue("grain_data_directory", std::string ()))[0];

  // shortcuts
  std::vector<boost::shared_ptr<emitting_grain> >& grains = 
    grain_collection_grain_model<chromatic_policy, T_rng_policy>::grains_;
  std::vector<array_1>& dnda =
    grain_collection_grain_model<chromatic_policy, T_rng_policy>::dnda_;

  // copy units from those supplied into *scatterer* unit map.
  T_unit_map& u = 
    mcrx::scatterer<chromatic_policy, T_rng_policy>::units_;
  u["length"] = units.get("length");
  u["wavelength"] = units.get("wavelength");
  u["mass"] = units.get("mass");

  // The grains use their preferred units, which are also copied into
  // the grain model units.

  // the grains are resampled onto a common grid, which is the first
  // grain's wavelengths. Since the PAHs have higher-resolution grids,
  // we put those first.
  const std::string pahnfile(p.defined("use_dl07_opacities") && 
			     p.getValue("use_dl07_opacities", bool()) ? 
			     "/PAH_neutral_DL07.fits" : "/PAH_neutral.fits");
  const std::string pahifile(p.defined("use_dl07_opacities") && 
			     p.getValue("use_dl07_opacities", bool()) ? 
			     "/PAH_ionized_DL07.fits" : "/PAH_ionized.fits");
  grains.push_back
    (boost::shared_ptr<emitting_grain>
     (new thermal_equilibrium_grain
      (data_directory + pahnfile, p)));
  grains.push_back
    (boost::shared_ptr<emitting_grain>
     (new thermal_equilibrium_grain
      (data_directory + pahifile, p)));
  // we must make sure we cut off the graphite grains at the largest
  // size of the pah's
  grains.push_back
    (boost::shared_ptr<emitting_grain>
     (new thermal_equilibrium_grain
      (data_directory + "/graphite.fits", p,
       grains[0]->sizes().back()*1.00001)));
  grains.push_back
    (boost::shared_ptr<emitting_grain>
     (new thermal_equilibrium_grain
      (data_directory + "/suv_silicate.fits", p)));

  T_unit_map& u2 = 
    mcrx::grain_model<chromatic_policy, T_rng_policy>::units_;
  u2["length"]=grains.front()->units().get("length");
  u2["wavelength"]=grains.front()->units().get("wavelength");
  // grains have no mass unit. But the size distributions do, so for
  // simplicity we select the scatterer mass unit as our mass unit.
  // NOTE that this arrangement means the size distributions and the
  // grains in general use different length units.
  u2["mass"] = u.get("mass");

  //  size distributions are made with our set *scatterer* units
  const std::string wd01_set =
    p.getValue("wd01_parameter_set", std::string());
  wd01_graphite_distribution graphite_size_distribution(wd01_set, units);
  wd01_silicate_distribution silicate_size_distribution(wd01_set, units);

  dnda.resize(grains.size());
  dnda[0].resize(grains[0]->asizes().size());
  dnda[1].resize(grains[1]->asizes().size());
  dnda[2].resize(grains[2]->asizes().size());
  dnda[3].resize(grains[3]->asizes().size());

  // evaluate the size distributions. Note that we need to convert
  // from grain length unit here to the scatterer units
  const T_float gsize_conv = units::convert(u2.get("length"),
					    u.get("length"));
  // ionized fraction is a simple fit to Fig 8 in DL07.
  const array_1 x_ion(0.55*log10(units::convert(u2.get("length"),"angstrom")
				 *grains[0]->asizes())-0.1);
  ASSERT_ALL(x_ion>=0);
  ASSERT_ALL(x_ion<=1);
  dnda[0] = (1-x_ion)*
    graphite_size_distribution.dn_da(grains[0]->asizes()*gsize_conv);
  dnda[1] = x_ion*
    graphite_size_distribution.dn_da(grains[1]->asizes()*gsize_conv);
  dnda[2] = graphite_size_distribution.dn_da(grains[2]->asizes()*gsize_conv);
  dnda[3] = silicate_size_distribution.dn_da(grains[3]->asizes()*gsize_conv);
  ASSERT_ALL(dnda[0]==dnda[0]);
  ASSERT_ALL(dnda[1]==dnda[1]);
  ASSERT_ALL(dnda[2]==dnda[2]);
  ASSERT_ALL(dnda[3]==dnda[3]);
}


template <template<class> class chromatic_policy, typename T_rng_policy> 
class mcrx::wd01_grain_creator_functor
{
public:
  typedef boost::shared_ptr<mcrx::grain_model<chromatic_policy, T_rng_policy> > return_type;
  
  wd01_grain_creator_functor() {};
  return_type operator()(const Preferences& p) const {

    //THIS IS A HACK BECAUSE WE CAN ONLY HAVE ONE ARGUMENT
    T_unit_map units;
    units["mass"] = p.getValue("massunit", std::string());
    units["length"] = p.getValue("lengthunit", std::string());
    units["wavelength"] = p.getValue("wavelengthunit", std::string());

    return return_type
      ( new mcrx::wd01_grain_model<chromatic_policy, T_rng_policy> 
	(p, units));
  };
};


#endif
